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Approach for Input Uncertainty 
Propagation and Robust Design 
in CFD Using Sensitivity 
Derivatives r 

An implementation of the approximate statistical moment method for uncertainty propa- 
gation and robust optimization for quasi I-D Euler CFD code is presented. Given uncer- 
tainties in statistically independent . random , normally distributed input variables , first - 
and second-order statistical moment procedures are performed to approximate the 
uncertainty in the CFD output. Efficient calculation of both first- and second-order sen- 
sitivity derivatives is required. In order to assess the validity of the approximations, these 
moments are compared with statistical moments generated through Monte Carlo simula- 
tions. The uncertainties in the CFD input variables are also incorporated into a robust 
optimization procedure. For this optimization, statistical moments involving first-order 
sensitivity derivatives appear in the objective function and system constraints. Second- 
order sensitivity derivatives are used in a gradient-based search to successfully execute a 
robust optimization. The approximate methods used throughout the analyses are found to 
be valid when considering robustness about input parameter mean values. 
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1 Introduction 

Gradient-based optimization of complex aerodynamic configu- 
rations and their components, utilizing high-fidelity Computa- 
tional Fluid Dynamics (CFD) tools, continues as a very active 
area of research (see, for example, [1-3]). In most of the CFD- 
based aerodynamic optimization and design studies to date, the 
input data and parameters have been assumed precisely known; 
we refer to this as deterministic or conventional optimization. 
When statistical uncertainties exist in the input data or parameters, 
however, these uncertainties affect the design and therefore must 
be accounted for in the optimization. Such optimizations under 
uncertainty have been studied and used in structural design disci- 
plines (see, for example [4-8]); we refer to these as nondetermin- 
istic or robust design optimization procedures. 

Derivatives of code output with respect to code input and pa- 
rameters are called sensitivity derivatives and they contain infor- 
mation which can be used to direct an optimization search. In a 
fluid flow optimization problem, the objective and constraint gra- 
dients are functions of the CFD sensitivity derivatives. Such de- 
rivatives can also be used to accurately approximate the CFD 
output in a small region, such as that near the mean value of a 
random variable. In [9], it is shown that a statistical First Order 
Second Moment (FOSM) method and automatic differentiation 
can be used to efficiently propagate input uncertainties through 
finite element analyses to approximate output uncertainty. This 
uncertainty propagation method is demonstrated herein for CFD 
code. 

An integrated strategy for mitigating the effect of uncertainty in 
simulation-based design is presented in [10]; this strategy consists 
of uncertainty quantification, uncertainty propagation, and robust 
design tasks or modules. Two approaches are discussed there for 
propagating uncertainty through sequential analysis codes: an ex- 
treme condition approach and a statistical approach. Both ap- 
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proaches can be efficiently implemented using sensitivity deriva- 
tives. For CFD code, the former approach is demonstrated in [ 1 1 ]. 
whereas the latter approach is demonstrated herein using second 
moment approximations. These uncertainty propagation methods 
have been developed and are being investigated as an alternative 
to propagation by direct Monte Carlo simulation lor potentially 
expensive CFD analyses. 

The present paper shows how the approximate statistical second 
moment methods. FOSM and the Second Order Second Moment 
(SOSMI counterpart, can be used in conjunction with sensitivity 
derivatives to propagate input data uncertainties through CFD 
code to estimate output uncertainties. The FOSM approximation is 
then used to perform sample robust optimizations. For demonstra- 
tion purposes, the method is illustrated on a simple design ex- 
ample which contains the important elements of more complex 
design problems. We assume that the input uncertainty quantifica- 
tion is given by independent norma Uy-distributed random vari- 
ables. and we demonstrate the strategy of [10] as applied to a CFD 
code module. This strategy is also applicable to correlated and/or 
non-normal ly distributed variables; however, the analysis and re- 
sulting equations become much more complex. 

The gradient-based robust optimization demonstrated herein re- 
quires second-order sensitivity derivatives from the CFD code. In 
[12], we present, discuss, and demonstrate the efficient calculation 
of second-order derivatives from CFD axle using a method pro- 
posed. but not demonstrated, in [13]. This method, used herein, 
incorporates first-order derivatives obtained by both forward- 
mode and reverse-mode differentiation in a noniterative scheme to 
obtain second-order sensitivity derivatives. 

To date, the only other demonstration or application of 
gradient-based, robust optimization involving advanced or high- 
fidelity (nonlinear! CFD code that we have found was just re- 
cently presented in [14-15]. The analytical statistical approxima- 
tion of their objective function for robust optimization also 
required second-order sensitivity derivatives. However, these 
studies employed a direct numerical random sampling technique 
to compute expected values at each optimization Mep in order to 


avoid the second derivatives. An example of linear aerodynamics 
involved in multidisciplinary performance optimization subject to 
yncertainty is found in [16]. 

XVo other aspects need to be pointed out in regard to the robust 
optimization demonstrations for CFD code modules presented 
herein and also in [14,15]. First, the sources of uncertainty con- 
sidered were only those due to code input parameters involving 
geometry and/or flow conditions; i.e., due to sources external to 
The CFD code simulation. Other computational simulation uncer- 
tainties, such as those due to physical, mathematical, and numeri- 
cal modeling approximations [17,18] — essentially internal model 
error and uncertainty sources, were not considered. That is, the 
discrete CFD code analysis results were taken to be deterministi- 
cally “certain” herein. Ultimately, all of these modeling sources 
of error and uncertainty must be assessed and considered. Sensi- 
tivity derivatives can also aid in this assessment [19] since the 
adequacy of an internal model's (i.e., algorithm, turbulence, etc.) 
prediction capability generally depends, to some extent, on the 
modeling parameter values specified as input. 

Second, as discussed in [15], uncertainty classification with re- 
spect to an event’s impact (from performance loss to catastrophic) 
and frequency (from everyday fluctuation to extremely rare) sets 
the problem formulation and solution procedure. Structural reli- 
ability techniques typically deal with risk assessment of infrequent 
but catastrophic failure modes, identifying the most probable 
point of failure and its safety index. Here, we are addressing the 
assessment of everyday operational fluctuations on performance 
loss, not catastrophe. Consequently, we are most concerned with 
aero performance behavior due to probable fluctuations. i.e.. near 
the mean of probability density functions (pdf). Structural reliabil- 
ity assessment is most concerned with improbable catastrophic 
events, i.e., probability in the tails of the pdf. Simultaneous con- 
sideration of both types of uncertainty is discussed in [16]. 

2 Integrated Statistical Approach 

In [10] an integrated methodology for dealing with uncertainty 
in a simulation-based design is proposed and demonstrated for a 
linkage mechanism design. The integrated strategy of [10] for 
mitigating the effect of uncertainty includes (a) uncertainty quan- 
tification. (b) uncertainty propagation, and (c) robust design. The 
present study utilizes the strategy proposed in [10], but differs in 
regard to uncertainty propagation and application. Here, we are 
considering the influence of uncertainty in CFD code input: that 
is. the effect of uncertainty in input geometry on aerodynamic 
shape-design optimization and the effect of uncertainty in flow 
conditions on design for flow control. 

2.1 Uncertainty Quantification. In this study, we consider 
the influence of uncertainty in CFD input parameterization vari- 
ables. We have assumed that these input variables are statistically 
independent, random, and normally distributed about a mean 
value. This assumption not only simplifies the resulting algebra 
and equations, but also serves to quantify input uncertainties. Fur- 
thermore, it is not an unreasonable assumption for input geometric 
variables subject to random manufacturing errors nor for input 
flow conditions subject to random fluctuations. 

2.2 Uncertainty Propagation. Uncertainty propagation is 
accomplished by approximate statistical second moment methods 
[9] and [20] where the required sensitivity derivatives are ob- 
tained by hand or by automatic differentiation. The first step in 
second moment analyses is to approximate the CFD system output 
solutions of interest in Taylor series form, truncated to the desired 
order. These approximations are formed to estimate the output 
value for small deviations of the input. 

Given input random variables b={£, b H } with mean 

values, b={£, B n }, and standard deviations. 

<r h J. first- and second-order Taylor series ap- 
proximations of the CFD output function. F. are given by 


and 
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respectively. Both first and second derivatives are evaluated at the 
mean values, b. 

One then obtains expected values for the mean (first moment) 
and variance (second moment) of the output function, F, which 
depend on the sensitivity derivatives and input variances, ar \ . The 
mean of the output function, F, and standard deviation <r F , are 
approximated as 
Rrst Order: 


F=F(b) 




(3) 


Second Order: 
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where both first and second derivatives arc evaluated at the mean 
values, b. Note in Eq. (4) that the second-order mean output. F. is 
not at the mean values of input b. i.e.. F^F(b). 

Equations (3) and (4) are the FOSM and SOSM approxima- 
tions, respectively, for the uncertainty propagation. The methods 
are straightforward with the difficulty largely lying in computation 
of the sensitivity derivatives. The very efficient and effective 
method used here to obtain such derivatives is presented in [12]. 


23 Robust Design. Conventional optimization for an ob- 
jective function. Obj, that is a function of the CFD output. F, stale 
variables. Q. and input variables, b. is expressed in Eq. (5). 
Herein, the CFD state equation residuals. R. are represented as an 
equality constraint, and other system constraints, g. are repre- 
sented as inequality constraints. The input variables, b. are pre- 
cisely known, and all functions of b are therefore deterministic, 
min Obj. 

Obj-Obj( F.Q.b) 

subject to 

R(Q,b) = 0 

g( F.Q.b 1=^0 (5) 

For robust design, the conventional optimization, Eq. (5), must 
be treated in a probabilistic manner. Given uncertainty in the 
input variables, b. all functions in Eq. (5) are no longer deter- 
ministic. The design variables are now the mean values, 
b={5j 5,,}, where all elements of b are assumed statisti- 

cally independent and normally distributed with standard devia- 
tions ctj, . The objective function is cast in terms of expected val- 
ues and becomes a function of F and The state equation 
residual equality constraint. R. is deemed to be satisfied at the 
expected values of Q and b. that is the mean values Q and b for 
the FO approximation. The other constraints are cast into a proba- 
bilistic statement: the probability that the constraints are satisfied 




is greater than or equal to a desired or specified probability, P k . 
This probability statement is transformed [10] into a constraint 
involving mean values and standard deviations under die assump- 
tion that variables involved are normally distributed The robust 
optimization can be expressed as 
min Obj, 

Obj = Obj(F, a F , Q,b) 

subject to 

R(Q,b)=0 

g(F,Q,b)+*<v*0, (6) 

where k is the number of standard deviations, <r g , that the con- 
straint g must be displaced in order to achieve the desired or 
specified probability, P k . For the FOSM approximation, standard 
deviations cr ¥ and sr % are of the form given in Eq. (3) involving 
first-order sensitivity derivatives. Therefore, a gradient-based op- 
timization will then require second-order sensitivity derivatives to 
compute the objective and constraint gradients. Note that for the 
SOSM approximation, third-order sensitivity derivatives would be 
required for these gradients. 

The calculation of second-order sensitivity derivatives for CFD 
code, such as those requited for SOSM and robust optimization 
with FOSM, was demonstrated in [13]; the efficient calculation 
method used herein is demonstrated and discussed in [12]. Both 
hand differentiation and automatic differentiation via the ADIFOR 
tool [21-23] were used in [12] and [13]. Both conventional and 
robust optimizations were performed using the Sequential Qua- 
dratic Programming method option in the Design Optimization 
Tools, DOT [24]. 

3 Application to Quasi 1-D Euler CFD 

A very simple example has been chosen to demonstrate the 
propagation of input uncertainty through CFD code and its effect 
on optimization. The quasi one-dimensional Euler equations and 
approximate boundary conditions are solved in a discretized CFD 
form (see Appendix). Of course, modem CFD methods are not 
required to solve the quasi l-D Euler equations; exact analytical 
solutions are available in any basic textbook on gas dynamics. 
However, these quasi I -D equations are commonly and effectively 
used as the initial test platform in the development of CFD algo- 
rithms and methods. 

Two separate applications are presented; the first involving 
propagation of geometric uncertainties, the second involving 
propagation of flow parameter uncertainties. Both uncertainty 
analyses are performed with quasi one-dimensional CFD Euler 
equations and boundary conditions describing subsonic flow 
through a variable area nozzle. The nozzle inlet is located at 
x~0 with area A(jt = 0)= 1 ; the nozzle outlet is at .v = 1. The area 
distribution is given by 

A(.v)= 1 - ax+bx 2 . 

The volume. V, occupied by the nozzle, is the integration of A(.r) 
over the length x—0 to x= 1 


obtaining and using the sensitivity derivatives needed for 
gradient-based optimization [25,26]. Therefore, we chose to by- 
pass issues related to this supersonic flow non-smoothness in 
these initial demonstrations of the present approach for uncer- 
tainty propagation and robust design for CFD code modules. 

3.1 Geometric Uncertainty Propagation. For the discus- 
sion of geometric uncertainty propagation, geometric shape pa- 
rameters a and b will represent the statistically independent ran- 
dom input variables, b. Die Mach number distribution through the 
nozzle, M, is viewed here as a component of the state variable, Q; 
its value at the inlet, M, is the CFD output, F. Applying the 
approach previously outlined (recall Eqs. (3) and (4)) yields the 
following first-and second-order approximations of the output 
function, M. 

Input random variables: b={a,6} 

CFD output function: F={M} 

First-order Taylor series: 

_ m m 

M(a,6) = M(a,6) + — (a~a)+ — (b-b) (7) 

da do 

Second-order Taylor series: 

_ dM _ d 1 M 

M (a,6) = Af(a,6) + — (a — a)+ -rr{b~b)+ ~r — (a- a) 
da do aadb 

X(Z7-5) + 0.5(^(«-t7)-|+0.5(^(fc-/7) 2 | 

( 8 ) 

The mean. M, and standard deviation a M of the output function 
are expressed as 
FOSM: 


M = M(</,6) 

(9) 



Predictions of M {a.b). M and <r M for FOSM (Eqs. (7) and (9)) 
and SOSM (Eqs. (8) and 00)) are compared with CFD solutions 
and Monte Carlo analyses based on CFD solutions, as given and 
discussed in the results section. 

3.2 Robust Shape Optimization. Applying the conven- 
tional optimization previously described yields 
min Obj. 


where a and b are the input geometric parameters. Three flow 
parameters are specified as input boundary conditions: the stagna- 
tion enthalpy, inlet entropy, and outlet static (back) pressure. The 
discretized quasi I-D Euler equation set is symbolically written as 
the state equation in Eq. (5); its residual. R is driven to (machine) 
zero for a solution (see Appendix). 

For supersonic flow through a variable area nozzle, shock 
waves generally appear and the flow solution (objective, con- 
straint. etc.) becomes noisy or non-smooth (see [25] ami the ref- 
erences cited therein). Care must be exercised with respect to 


Obj=Obj(M M.b) 

subject to 

R(M.«.b)-0 

V(a.fr)<a 


(ID 


where the system constraint. V . is a constraint on the nozzle vol- 
ume and depends only on a and b: our objective does not explic- 
itly depend on M. 

Applying the robust optimization previously described yields 
min Obj. Obj= Obj( filer v .<7.6 ) 

subject to 
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R ( M , ff ,£)~0 

V(a$)+ka v ^ 0 , 



With a and b subject to statistical uncertainties (which may be 
due to measurement, manufacturing, etc.), V becomes uncertain. 
Since V is linearly dependent on a and b , it is also normally 
distributed. Therefore, its standard deviation, <t v , is given exactly 

by Eq. (13)- 

To demonstrate the optimizations, a simple target-matching 
problem is selected; a unique answer is obtained when an equality 
volume constraint is enforced. The CFD code is run for given a 
and b\ the resulting M (a,b) and corresponding V(a,b) are taken 
as the target values Mt and Vt, respectively. For this conventional 
optimization, the objective function and constraint function for V 
of Eq. (11) become 

Obj(M,a,£) = [M(a,&)~ Mt] 2 

V(a,b)-\ t=0 

enforced as 


V(a,b)- Vt^O and Vx-V(a,b)**0 (14) 

for the convenience of the optimizer 
For robust optimization using the FOSM approximation, the 
corresponding objective and constraint on V of Eq. (12) become 

Obj( M, t r v ,a,b) = [M(a,/>)- Mt] 2 + cr 2 ^ 

V(aj)-Vt+k<r v =0 

similarly enforced as 

V(a,b)-Vt+ka v ^Q and Vt- V(a,b)-k<r v **Q. (15) 

Note that for <T a =<r h = 0 in Eq. (15), the conventional optimiza- 
tion is obtained. Also, in the probabilistic statement of the con- 
straint on V, it is assumed that the desired volume is less than or 
equal to Vt. 

33 Flow Parameter Uncertainty Propagation. A second 
example of uncertainty in CFD involves fluctuations in input flow 
parameters. For the discussion of flow parameter uncertainty 
propagation, the free-stream Mach number, Minf. and the nozzle 
static back pressure. Pb. will be taken as statistically independent 
random variables. Specifying the free-stream Mach number sets 
the stagnation enthalpy. The Mach number distribution through 
the nozzle. M, is again v iewed as a component of the state vari- 
able. Q; its value at the inlet, M, is the CFD output F. Applying 
the approach previously outlined yields equations which are simi- 
lar to Eq. (7) through (10) but with 

Input random variables: b= {Minf.Pb} 

CFD output function: F={M} 

Again, predictions of M, M, and <r M for the FOSM and SOSM 
approximations are compared with CFD solutions and Monte 
Carlo analyses based on CFD solutions, as given and discussed in 
the next section. 

3.4 Robust Design for Flow Control. The conventional 
optimization is expressed as 

min Obj Obj- Obj(M. Minf.Pb) 

subject to 


R(M,Minf.Pb)=0 

<?(Minf,Pb)<0. 


U6) 


w here q is a constraint on the mass flux through the nozzle. 


The robust optimization is expressed as 

min Obj, Obj=Obj(M,o’ W ,Minf,Pb) 
subject to 

R( M, Minf, Pb) = 0 

(17) 

q ( Minf, Pb) + k <r q ^ 0. 

For the free-stream Mach number, Minf, and the nozzle back 
pressure, Pb, subject to statistical uncertainties, the mass flux, q , 
becomes uncertain. Since q is dependent on Minf and Pb, its stan- 
dard deviation, a q , may be approximated by 



Since q is not a linear function of Minf and Pb, Eq. (18) is not 
exact 

To demonstrate the optimizations, a simple target-matching 
problem is again chosen. The CFD code is run for given Minf and 
Pb, the resulting M and corresponding q are taken as the target 
values Mt and qt , respectively. For this conventional optimization, 
the objective function and constraint functions of Eq. (16) are 

Obj(M.Minf,Pb) = [M(Minf,Pb) - Mt] 2 

< 7 (Minf,Pb)-<// = 0 

enforced as 

qi Minf.Pb) — and <// -<y( Minf.Pb) ^0. (19) 

For robust optimization using the FOSM approximation, the cor- 
responding objective and constraint on q of Eq. (17) can be shown 
as 


Obj=Obj(M.<r A/ .Minf,Pb) = [M(Minf,Pb)-Mt] 2 + o- 2 f 
q{ Minf.Pb ) — qt+k tr q - 0 

enforced as 

q(Minf.Pb)-f//+4*<r^() and qt -q( Minf.Pb) -ktr^Q. 

(20) 

Again note that for fr NUlll =(r ln> =0 in Eq. (20), the conventional 
optimization is obtained. Also, in the probabilistic statement of the 
constraint on q , it is assumed that the desired mass flux is less than 
or equal to qt. 


4 Sample Results and Discussion 

Presentation and discussion of results for the sample quasi l-D 
Euler CFD problems are divided into four topics: function ap- 
proximations. uncertainty propagation, pdf approximations, and 
robust optimization. For the first three topics, the FOSM and 
SOSM approximations are assessed by comparison with direct 
CFD simulations. Note that the CFD code and its corresponding 
derivative codes are executed at the mean values of the random 
input variables to obtain the functions and sensitivity derivatives 
needed for construction of the statistical approximations. 

4.1 Function Approximations. It is important to assess the 
Taylor series output function approximations with direct nonlinear 
CFD code simulations prior to presenting uncertainty propagation. 
If the CFD output function. M. is quasi-linear with respect to the 
input variables of interest, one can expect first-order approxima- 
tions to be reasonably good; that is. the first-order moments given 
by Eq. (3) should match well with the moments produced by a 
Monte Carlo simulation. For a more nonlinear system, one natu- 
rally expects better accuracy with second-order approximations; 
that is. uncertainty analyses which include second-order terms 
should yield results which better predict the statistical moments 
produced by the Monte Carlo simulation. 
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Fig. 1 Comparison of function approximations versus CFD 
solution, input variable b=b 


Figures 1-4 show that for F=M(c/,/>), M behaves as a quasi- 
linear function in the neighborhood of (n,6), whereas for F 
= M(Minf.Pb), M is more nonlinear in the neighborhood of 
(Minf.Pb). In these figures, approximations of the CFD output 
functions. M (a,b) and M(Minf.Pb), using the first- (FO) and 
second-order (SO) Taylor series (as given in Eqs. (7) and (8) for 
M(a,6)), are compared to direct solution of the Euler CFD. In 
each example, two traces were considered through the design 
space. Trace 1 varied the first input variable, while the second 
remained fixed at its mean value, and vice versa for trace 2. The 
required first- and second-order sensitivity derivatives needed for 
construction of the first-and second-order approximations were 
obtained by hand differentiation and automatic differentiation as 
discussed and presented in [12]. 

Nonlinear behavior of the CFD result is reasonably well ap- 
proximated by the SO result in all plots; however, there does 
appear to be an inflection point in the CFD results given in Fig. 3. 
Note that the linear FO result is a good approximation in the 



Fig. 2 Comparison of function approximations versus CFD 
solution, Input variable 
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Fig. 3 Comparison of function approximations versus CFD 
solution, input variable Pb=Pb 


geometric example; the flow parameter example is more nonlin- 
ear. At larger deviations from the mean, a linear approximation for 
M(Minf.Pb) loses accuracy. 

4.2 Uncertainty Propagation. Approximation of the statis- 
tical first and second moments is done using equations £q. (9) and 
(10) (geometric example), and corresponding equations for the 
flow parameter example. Again, both first- and second-order sen- 
sitivity derivatives are required and the prediction is straightfor- 
ward. given these derivatives. An independent verification of these 
approximate mean and standard deviation values is obtained here 
using direct Monte Carlo (MC) simulation with the quasi l-D 
Euler CFD code and standard statistical analyses of these Monte 
Carlo results. The standard statistical analyses used were from 
Microsoft <S> Excel 20(X) and the random number generator 
MZRAN used was from [27]. Tables I and 2 give results for the 
mean (first moment) and Tables 3 and 4 give results for the stan- 



Fig. 4 Comparison of function approximations versus CFD 
solution, input variable Minf-Minf 
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Ifcble 1 Percent difference from MC for FO and SO predictions 
of A (1,6) 


Case 

Input <r 

<V=°b 

M 

MC 

% Error 
MC 

% diff w/MC 
FO Predict 

% diff w/MC 
SO Predict 

' 1 

0.01 

0.4041 

0.0187 

-0.0105 

0.0656 

2 

0.02 

0.4040 

0.0379 

0.0716 

0.1531 

3 

0.04 

0.4054 

0.0756 

-0.2867 

0.0383 

4 

0.06 

0.4055 

0.1142 

-0.3012 

0.4301 

5 

0.08 

0.40% 

0.1557 

-1.3078 

-0.0209 


Table 2 Percent difference from MC for FO and SO predictions 
of M (Minf.Pb) 

Case 

Input a 

^MwT^Pb 

M 

MC 

% Error 
MC 

* diff w/MC 
FO Predict 

% diff w/MC 
SO Predict 

1 

0.0 1 

0.3933 

0.0056 

0.0037 

-0.0269 

2 

0.02 

0.3932 

0.0114 

0.0187 

-0.1034 

3 

0.04 

0.3898 

0.0229 

0.8917 

0.3991 

4 

0.06 

03889 

0.0364 

1.1251 

0.0141 


Table 3 

Of <7||, 

Percent difference from MC for FO and SO predictions 
geometric example 


Input a 

O’M 

% diff w/MC 

% diff w/MC 

Case 

* 

H 

MC 

FO Predict 

SO Predict 

i 

0.0! 

0.0102 

-0.5773 

-0.5708 

2 

0.02 

0.0207 

-1.7026 

-1.6769 

3 

0.04 

0.0414 

-1.5794 

-1.4766 

4 

0.06 

0.0625 

-2.2590 

-2.0296 

5 

0.08 

0.0853 

-4.3987 

-4.(XX)I 


Table 4 Percent difference from MC for FO and SO predictions 
of <r M , flow parameter example 

Case 

Input <r 

a Mmf~ ^Ph 

MC 

diff w/MC 
FO Predict 

% diff w/MC 
SO Predict 

1 

0.0 1 

0.0030 

1.1815 

1.2473 

2 

0.02 

0.0062 

-1.5093 

-1.2533 

3 

0.04 

0.0125 

-4.1604 

-3.1680 

4 

0.06 

0.0199 

-4.4070 

-2.1938 


dsri deviation (second moment) value comparisons. The input 
deviations (<r a and a b ) or ((Ty^s and o> 6 ), are taken to be equal 
md given in the second column of each table. The third column in 
each table gives the result from the MC simulation, where the 
sample size (N) used was 3000. The MC error in its predicted 
mean is <r M / JV. which is given in the fourth column of Tables 1 
and 2. The first-order (FO) and second-order (SO) approximate 
predictions are given in the last two columns of each table as 
percent difference from the MC results. 

Results from the last two columns of Tables 1-4 are plotted in 
Figs. 5 and 6. Note that mean and standard deviation approxima- 
tions for output M agree well with the MC results for small stan- 
dard deviations in input variables, (<r a ,<r b ) or for 

both the FO and SO predictions. At higher standard deviations, the 
SO results tend to agree better with the MC results, especially in 
the flow parameter example. 

As expected, the SO approximations are generally better than 
those for FO and an increase in the standard deviation of input 
parameters monotonically increases the error associated with both 
predictions. Both FO and SO approximations predict the mean 
more accurately than the standard derivation. Also note that the 
MC simulation with a sample size of 3000 limits accuracy of the 
MC results. This is apparent in the raggedness of the pdfs shown 
in the following section. 

43 Probability Density Function Approximations. Given 
a mean and standard deviation of the CFD output function (from 
either an MC simulation or an FO or SO prediction) and assuming 
a normal distribution, one may then construct a pdf to approxi- 
mate the behavior of the nondcterministic output function. This 
approximation is compared to pdf histograms generated from MC 
simulations in Figs. 7. 8. and 9. The bars depict the MC histo- 
gram. and the solid curve is a normal distribution at the MC mean 
value and MC standard deviation as given in the previous tables. 
The MC simulation size of 3(XX) is certainly not sufficient to ob- 
tain a smooth pdf. We note that both the first-order and second- 
order normal distributions are indistinguishable from this normal 
MC curve at this scale so neither is shown in the figure. !t is 
^parent that for either the quasi-iinear functional dependence on 
a and b (Fig. 7). or for small input standard deviations of the flow 
parameters (Fig. X). the statistical approximations arc good for a 
significant region about the mean but tend to break down in pre- 
dicting the tails of the distribution. This is significant, for if one is 
primarily interested in reliable failure predictions, us for structural 
design, this prediction may not be good enough. It is felt, how- 



Fig. 5 Comparison of statistical moment approximations with Monte Carlo 
simulation results, geometric examples 




Fig. 6 Comparison of statistical moment approximations with Monte 
Carlo simulation results, flow parameter examples 


ever, that in aerodynamic performance optimization using CFD, 
where robustness about the mean is desired, these approximations 
may be good enough. 

It is not surprising that a nonl inear CFD output function be- 
haves differently for randomness in different input variables. For 
M(Minf,Pb) at higher input standard deviations (cr Millf =a> fr 


(Input stdvs=O06 output sttfvsa06S30) 
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Fig. 7 Probability density function for M(a,t>) for <r M ^(r b 
-0.08 


(input stdvs^LOS output std«4tj012Z78) 

aoe-i — 
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=0.06), the pdf of the output function is no longer normal. In Fig. 
9 one can see the non-normal behavior of CFD output given nor- 
mally distributed input variables Minf and Pb. 

4.4 Robust Optimizations. Optimization results were gen- 
erated using the quasi I-D Euler CFD code and the procedure 
given by Eqs. (5) and (6). As noted earlier, conventional optimi- 
zation is obtained for (r a =a h = 0 or for <r Min r=fr Ph =0. For the 
FOSM approximation, first-order sensitivity derivatives are re- 
quired to obtain (r Xf and tr v (Eqs. (9) and (13)) or & q (Eqs. (9) and 
(18)); therefore, second-order sensitivity derivatives will be re- 
quired for their gradients, which are used in the optimization. The 
second-order derivatives were obtained by the manner presented 
in [12], except for those associated with Eq. (18) where the first- 
order derivatives were finite differenced simply for convenience. 

It is seen from Eqs. (12) and (17) that the robust optimization 
results should depend on the probabilistic parameters Ur a% (r h ) or 
(°Minf ,o>,,). and k. The desired probability, P k , is that from the 
normal cumulative distribution function since <r v and here are 
assumed to be normally distributed. For each robust optimization 
example, two cases are presented. For case I, P k is fixed at 
k~ I, i.e.. Pi =84.13#, and the effect of increasing the in- 
put variable standard deviations is addressed. For case 2, the stan- 
dard deviations of the input variables are fixed at 0.0! and P k 
increases. 



Fig. 8 Probability density function for M (Mlnf,Pb) for Fig. 9 Probability density function for M (Minf,Pb) for 

**<rpfc"0.Q2 *<Tp 6 *0.06 
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Tab** 5 Robust shape optimization results with increasing in- We 6 Robust shape optimization results with increasing P k 
put parameter <r for k~1 for 


b 

I 

i 

Obj 

ft 


<r v 

K 

p* 

* 


Obj 

& 


<r v 

0.00 

0.6001 

0.3001 

0.0000 

0.4043 

0.0000 

0.0000 

0 

0.5000 

05996 

02995 

0.000104 

0.4041 

0.0101 

0.006 

0.02 

0.6685 

0.3667 

0.0004 

0.4036 

0.0203 

0.0120 

1 

0.8413 

0.6246 

0.3189 

0.000118 

0.4004 

0.0101 

0.006 

0.04 

0.7338 

0.4286 

0.0016 

0.4018 

0.0406 

0.0240 

2 

0.9772 

0.6698 

0.3687 

0.000104 

0.4041 

0.0101 

0.006 

0.06 

0.7948 

0.4841 

0.0037 

0.3984 

0.0607 

0.0360 

3 

0.9986 

0.7052 

0.4037 

0.000104 

0.4042 

0.0102 

0.006 

0.08 

0.8534 

0.5358 

0.0065 

0,3941 

0.0804 

0.0480 

4 

0.9999 

0.7406 

0.4388 

0.000104 

0.4043 

0.0102 

0.006 


4.4./ Robust Shape Optimization Results . In Table 5, results 
for case 1 of the robust shape optimization are displayed. For 
a a -cr b ranging from 0 to 0.08, optimal values for die input vari- 
ables (a,E) are listed. As a a —a b increases, so does o v . Accord- 
ingly. the mean values, (a.b). which minimize the objective func- 
tion and satisfy die probabilistic constraint, become increasingly 
displaced from the target volume, Vt. This is shown in Fig. 10. 
Mean values (a,B) change, keeping the mean value, M (u,£), of 
the probabilistic output near the taiget value, Mt. The robust de- 
sign points track the dashed curve for M=Mt with some displace- 
ment due to the <r M 2 term of the objective, Eq. U5). The volume, 
V(a.b), is displaced from the solid curve V=Vt by ka Vl as 
required by the probabilistic constraint. This displacement can be 
viewed as the solution dependent or “effective” safety factor. 

In Fig. 1 1 the changing area distribution of the robust optimi- 



Fig. 10 Optimization results in design space (a,b), P k fixed at 



zatton is illustrated. As the standard deviations of design variables 
(a,b) increase, the optimal nondeterministic volume, V(a,b), sig- 
nificantly decreases. 

The results for case 2 of the robust shape optimization, where 
a a =a b is fixed at 0.01, and P k increases from 50 percent to 99.99 
percent (k=0 to 4) are given in Table 6. Again mean values (a^b) 
change, keeping the mean value, M (a,b), of the probabilistic out- 
put near the taiget value, Mt. Since o a = cr b remains small, the o ^ 
term of the objective remains small, and the displacement of M 
from the dashed line depicting Mt due to the <r M term remains 
small as shown in Fig. 12. With an increase in P A , V(a^b) is 
displaced from the solid curve V- Vt by ko v , as required by the 
probabilistic constraint. Accordingly, the mean values, {ajb), 
which minimize the objective function and satisfy the constraint, 
again become increasingly displaced from those at the taiget vol- 
ume, Vt. Note the significant displacement of the solution from 
the taiget volume when P k is large, i.e., when one is attempting to 
incorporate the tails of the pdf. In order to increase the probability 
of constraint satisfaction from 97.77 percent to 99.99 percent, one 
sees a significant change in ( ajj ) for a mere gain of 2 percent in 
constraint satisfaction. 

4.4.2 Robust Design for Flow Control Results. Similar re- 
sults are seen in the flow parameter example. In Table 7, the 



Fig, 12 Optimization results in design space (a,b), <r fixed at 
0.01 


Table 7 Robust design for flow control results with increasing 
input parameter o* for k=1 


( r s* n —< T rb 

Minf 

Pb 

Obj 

M 

O'M 


0 

0.3000 

0.8000 

0.0000 

0.3933 

0.0000 

0.0000 

0.02 

0.2861 

0.7883 

0.0001 

0.3974 

0.0116 

0.0058 

0,04 

0.2655 

0.7801 

0.0005 

0.3985 

0.0231 

0.0112 

0.06 

0.2555 

0.7653 

0.0012 

0.4050 

0.0327 

0.0163 

0.08 

0.2468 

0.7498 

0.0020 

0.4118 

0.0407 

0.0209 
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Fig. 13 Optimization results in design space (Minf,Pb), P k 
fixed at P, 


results for case 1 are displayed. For o- Minf =<7 Pb ranging from 0 to 
0.08, optimal values for the input variables (MinfjPb) are listed. 
As (r Mlnf = a Pb increases, so does a q . Accordingly, the mean val- 
ues, (Minf,Pb), which minimize the objective function and satisfy 
the constraint, become increasingly displaced from the target mass 
flux, qt. This is shown in Fig. 13. Mean values (MinfPb) change, 
keeping the mean value, M(Minf,Fb), of the probabilistic output 
near the target value, Mt. The robust design points track the 


Table 8 Robust design for flow control results with increasing 
P k for cr»0.01 


K 


Minf 


Obj 

M 

<*\l 

Oq 

0 

0.5000 

0.3000 

0.8000 

0.00003 

0.3933 

0.0060 

0.0030 

1 

0.8413 

0.2919 

0.7953 

0.00003 

0.3945 

0.0059 

0.0029 

2 

0.9772 

0.2825 

0.7916 

0.00003 

0.3949 

0.0059 

0.0029 

3 

0.9986 

0.2688 

0.7896 

0.00003 

03936 

0.0060 

0.0028 

4 

0.9999 

0.2598 

0.7867 

0.00003 

0.3938 

0.0060 

0.0028 



ftM-Sfrwm Uacfc No. (IM) 


Fig. 14 Optimization results in design space (Minf,Pb), a 
fixed at 0.01 


dashed curve for M=Mt with displacement due to the <7^ term of 
the objective, Eq. (17). The optimized mass flux, $(Minf,Pb), is 
displaced from the solid curve q=qt by kcr q as required by the 
probabilistic constraint 

The results for case 2 of the robust design for flow control, 
where (z Millf = a Pb is fixed at 0.01, and P k increases from 50 per- 
cent to 99.99 percent (*=G to 4) are given in Table 8. Again, 
mean values (Minf,Pb) change, keeping the mean value, 
M(Minf,Pb), of the probabilistic output near the target value, Mt. 
As in the preceding example, since o- Minf = remains small, the 
term of the objective remains small and the displacement due 
to the ar~ M term remains small, as shown in Fig. 14. With an 
increase in P k , </(Minf,Pb) is displaced from the solid curve q 
— qt by k<r q , as required by the probabilistic constraint Accord- 
ingly, the mean values, (Minf,Pb), which minimize the objective 
function and satisfy the constraint again become increasingly dis- 
placed from the taiget mass flux, qt . Again, note the significant 
displacement from the target mass flux incurred in the higher 
probability optimizations, i.e., when one is attempting to incorpo- 
rate the tails of the pdf. 

5 Concluding Remarks and Challenges 

The present results represent an implementation of the approxi- 
mate statistical moment method for uncertainty propagation and 
robust optimization for a quasi I-D Euler CFD code. Assuming 
statistically independent, random, normally distributed input vari- 
ables, first- and second-order statistical moment procedures were 
performed to approximate the uncertainty in the CFD output. Ef- 
ficient calculation of both first- and second-order sensitivity de- 
rivatives was employed and the validity of the approximations 
was assessed by comparison with statistical moments generated 
through Monte Carlo simulations. The uncertainties in the CFD 
input variables were incorporated into a robust optimization pro- 
cedure where statistical moments involving first-order sensitivity 
derivatives appeared in the objective function and system con- 
straints. Second-order sensitivity derivatives were used in a 
gradient-based robust optimization. The approximate methods 
used throughout the analyses were found to be valid when con- 
sidering robustness about input parameter mean values. 

Collectively, these results demonstrate the possibility for an ap- 
proach to treat input parameter uncertainty and its propagation in 
gradient-based design optimization that is governed by complex 
CFD analysis solutions. It has been demonstrated on a very simple 
CFD code and problem; there are computational resource issues to 
be addressed in application to significant 2-D and 3-D CFD codes 
and problems. Conventional optimization in 2-D and 3-D already 
requires significantly more resources than in 1-D and obtaining 
second-order sensitivity derivatives for more design variables will 
add even more work. Some of these are addressed in [12] and 
work is presently in progress regarding application and demon- 
stration using 2-D and 3-D Euler CFD code. 
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Nomenclature 

A = nozzle area 
a = geometric shape parameter 
b = geometric shape parameter 
b = vector of independent input variables 
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Transactions of the ASME 


F = vector of CFD output functions 
g = vector of conventional optimization constraints 
k = number of standard deviations 
M = Mach number at nozzle inlet 
M = vector of Mach number at each grid point 
Minf = free- stream Mach number 
Ml = target inlet Mach number 
N = Monte Carlo sample size 
Pb = normalized nozzle static back (outlet) pressure 
Q — vector of flow-field variables (state variables) 
q — mass flux through nozzle 
qt - target mass flux through nozzle 
R = vector of state equation residuals 
V — nozzle volume 

Vt = target nozzle volume used for optimization 
x = normalized axial position within nozzle 
cr = standard deviation 
a 2 = variance 


Superscript 

- = mean value 


Appendix-Quasl One-Dimensional Euler CFD 

The steady-state, quasi 1 -D Euler equations are 


mQ ) 

dx 




where 

Q=[p,pu,fte„] T , 

E(Q) = [pu.pu 2 + PApe n + P)u] T , 
dA I 

S(Q)=--jjj[pu,pu-Ape„+P)u] J . 


In these equations, p is density, it is flow speed. P is pressure, e„ 
is the specific total eneigy (i.e„ e {) — e + ir/2, where e is the spe- 
cific internal energy), and A(x) is the local cross-sectional area. 
The ideal gas law with a constant ratio of specific heats y (taken 
to be 1.4) is used for closure, which implies P~( y— I )(pe<, 
— p«'/2). The governing equations are discretized and solved nu- 
merically with the upwind flux-vector-splitting method of Van 
Leer [28]. which includes the use of higher order accuracy to 
approximate the flux terms. A more complete discussion of these 
numerical procedures is presented in [29]. The flow field was 
discretized with 100 grid points. This discretiztion of the govern- 
ing equations, together with the numerical treatment of the bound- 
ary conditions, results in a large set of coupled nonlinear algebraic 
equations with the form of Eq. (5). In this study, the procedure for 
solving the discrete nonlinear flow equations is Newton's method. 
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